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ABSTRACT 

The heating, acceleration, and pitch-angle scattering of charged particles by MHD turbulence are 
important in a wide range of astrophysical environments, including the solar wind, accreting black 
holes, and galaxy clusters. We simulate the interaction of high-gyrofrequency test particles with fully 
dynamical simulations of subsonic MHD turbulence, focusing on the parameter regime with /3 ~ 1, 
where j3 is the ratio of gas to magnetic pressure. We use the simulation results to calibrate analytical 
expressions for test particle velocity-space diffusion coefficients and provide simple fits that can be 
used in other work. 

The test particle velocity diffusion in our simulations is due to a combination of two processes: 
interactions between particles and magnetic compressions in the turbulence (as in linear transit- 
time damping; TTD) and what we refer to as Fermi Type-B (FTB) interactions, in which charged 
particles moving on field lines may be thought of as beads sliding along moving wires. We show that 
test particle heating rates are consistent with a TTD resonance which is broadened according to a 
decorrelation prescription that is Gaussian in time (but inconsistent with Lorentzian broadening due 
to an exponential decorrelation function, a prescription widely used in the literature) . TTD dominates 
the heating for v s ^> va (e.g. electrons), where v s is the thermal speed of species s and va is the 
Alfven speed, while FTB dominates for v s <C va (e.g. minor ions). Proton heating rates for j3 ~ 1 are 
comparable to the turbulent cascade rate. Finally, we show that velocity diffusion of collisionless, large 
gyrofrequency particles due to large-scale MHD turbulence does not produce a power-law distribution 
function. 



1. INTRODUCTION 

The interaction between charged particles and mag- 
netohydrodynamic (MHD) turbulence plays a role in 
the energy balance of diverse astrophysical environ- 
ments such as the solar corona and solar wind (e.g 
ICranmer et al.l 120071) and accretion d isks around black 
holes (e.g. iQuataert fc Gruzinovl 119991 ) . The coupling 
between turbulence and particles is also important for 
the t ransport and confinement of cosmic ray s in galaxies 
(e.g. !Chandradl2l)00t lYan fc Lazarian1[200l . 

This paper focuses on the interaction of test par- 
ticles with subsonic (and thus weakly compressible) 
MHD turbulence. Such turbulence consists primar- 
ily of nonlinearly inter acting Alfven waves which drive 
the turbulent cascade ( Kraichnan 1965; Shcbali n et al.l 
1 1983t IGoldreich fc Sridharll995l) . along with slow magne- 
tosonic m odes that are ad vected passiv el y by the Alfvenic 
casca de (jLithwick fc Goldreichl 120011: Cho fc Lazarianl 
120021 ). On observational (IChen et al.l I2011D. the- 



oretical fGoldreich fc Sridhar! 119951). and numerical 
(|Maron fc GoldreicHl2001l:lBeresnvakll20ill) grounds, the 
MHD cascade is believed to be strongly anisotropic, 
with most power in the inertial range of the cas- 
cade in modes with wavevectors primarily perpendic- 
ular to the magnetic field. This has many impli- 
cations for the coupling between particles and tur- 
bulence, including e.g. that cyclotron heating of 
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particles by the MHD c ascade is significantly sup 
press ed (IQuataertl 119981: ICranmer fc van Ballegooiie 7 
120031 t hough seelLeamon et al. l (ll999l ): IGarv fc Borovsk . 
(|2008l ): |jiang et al.l ([2009D for discussions of the cyclotron 
resonance in the corona and solar wind). 

The interaction between test particles and plasma 
waves has been ext ensively studied in the "quasilin- 
ear" approximation (jKennel fc Engelmandll9 66: Jokipii 
I1966D . which stipulates that test particles execute unper- 
turbed helical motion around magnetic field lines, and 
that plasma waves are long-lived relative to their pe- 
riods. This implies that wave-particle interactions and 
energ y exchange occur only at discrete resonances (IStixl 
1992, and references therein). To the extent that MHD 
turbulence is well-described by a superposition of long- 
lived small-amplitude plasma waves, quasilinear theory 
will accurately describe test particle diffusion and heat- 
ing in turbulence. 

However, the picture of strong anisotropic MHD tur- 
bulence de veloped over the last ~30 years by many au- 
thors (e.g. IMontgomerv fc Turner! |1981 | : Shcbali n et al.l 
Il98l IHigdonl 119841: IGoldreich fc Sridhar! 119951 hence- 
forth GS) suggests that the Alfven and slow waves com- 
prising weakly compressible MHD turbulence are not 
long-lived, and instead decorrelate due to non-linear in- 
teractions before they can propagate over distances of 
multiple wavelengthsQ As a result, the discrete reso- 
nances of quasilinear theory are expected t o be substan- 
tially broadened in MHD turbulence (e.g. iBieber et al.l 

4 The use of strong here refers to this state of "critical balance" 
between eddy and wave timescales, rather than the amplitude of 
the turbulence. 
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1991 IShalchi et all [200l IShalchi fe Schlickeiserl l200i 
Qin et al.ll2006t lYan fc Lazarianll2008D . 

In this paper, we study the interaction between test 
particles and driven MHD turbulence (see al so earlier 
work bv lDmitruk et al.ll2004t iLehe et al.ll2009f ). In par- 
ticular, we quantify the velocity-space diffusion, particle 
heating, and particle acceleration that results. We com- 
pare these numerical results in detail to analytic models 
and, in particular, calibrate models of resonance broad- 
ening. Our focus is on the basic physics — the diffusion 
coefficients we calculate (e.g. eqn. [10] & E] & Table [2]) 
can be used for a wide range of applications, some of 
which we will explore in detail in future studies. 

The paper is structured as follows. In <j2j we sum- 
marize the qualitative features of test particle interac- 
tions with turbulence, including both resonance broad- 
ening and non- resonant interactions. In $3] we use these 
physical ideas to derive analytical expressions for velocity 
diffusion coefficients in turbulence, while in 21 we calcu- 
late the resulting heating rates for a thermal distribution 
of test particles. Many of these results are not new, but 
they provide a useful analytic framework for interpreting 
our test particle numerical results and so are included for 
completeness. In <|5]we describe our numerical methods 
for evolving test particles in simulations of MHD turbu- 
lence, and in we compare our analytical predictions 
to the results of test particle simulations. Finally, in $7] 
we discuss the conclusions and implications of our work. 
We also include several Appendices which consider re- 
lated ideas. In Appendix El we discuss the interaction 
of test particles with one finite-amplitude wave, in Ap- 
pendix [B] we consider an extension of our model into a 
the regime of weak turbulence, and in Appendix [C] we 
discuss the power spectra of our turbulence simulations 
in the context of weak and strong turbulence. 

2. QUALITATIVE DISCUSSION OF TEST PARTICLE 
TRANSPORT IN MHD TURBULENCE 

In this paper, we focus on isothermal MHD turbulence 
with j3 = pCs/(Bq/8it) ~ 1 and e such that the turbu- 
lence is subsonic and sub-Alfvenic, where e is the en- 
ergy input rate per unit mass into the turbulence. Ad- 
ditionally, we focus on high-gyrofrequency particles with 
Q 3> Wmaxj where fl = qB/mc is the particle cyclotron 
frequency and w max is the maximum resolveable wave 
frequency in the turbulence. Our motivation for doing 
so is that this inequality is believed to be satisfied even 
deep in the inertial range of weakly co mpressible MHD 
turbulence (see e.g. iHowes et al.l l2008h . In addition to 

3> w max , magnetic moment conservation generally also 
requires that the amplitude of the turbulent fluctuations 
on scales of the Larmor radius r j, satisfy 5v/v\ <gC 1 
dMcChesnev et alJll987t lChastonll200l iChandran et al.1 
1201(11 . This is satisfied in our simulations both because 
we focus on subsonic f) ~ 1 turbulence and because the 
Kolmogorov power spectrum that the turbulence self- 
consistently develops has only modest power on small 
scales ~ r£. We shall see that magnetic moment con- 
servation is indeed reasonably well satisfied in our test 
particle simulations (see ^6|. This implies that particle 
acceleration will be primarily in the parallel direction. 
(We use the subscripts || and _L throughout to indicate 
parallel and perpendicular to the local magnetic field, 
respectively.) 



Generically, there will be two processes that cause 
changes in parallel velocity for the high-57 particles un- 
der consideration. The first is transit-time damping 
(TTD), which is analogous to Landau damping. In a 
spatially-varying magnetic field, charged particles feel 
mirror forces, given by /xVii-B, where /i = mv\j1B is 
the particle's magnetic moment. If the spatial variation 
is provided by a compressive wave moving with a phase 
speed v p , then particles with Uii ~ v p will "surf" the wave 
and experience correlated acceleration for long times (un- 
til the particle is accelerated such that it is no longer in 
resonance). On the other hand, a particle with vu v p 
will experience time- varying accelerations that will av- 
erage to zero over long times. Thus, over time scales 
sufficiently short that particle velocities do not change 
substantially (where linear theory applies) the interac- 
tion of the wave with a distribution of particles will be 
given by a delta-function, Dm cx 8(v\\ — Vp), where Dii is 
the parallel velocity diffusion coefficient. 

This picture will be modified for the interaction of par- 
ticles with strong turbulence. Strong turbulence can be 
thought of as a distribution of waves, which nonetheless 
do not propagate long distances as waves, but instead de- 
cohere on a timescale where u> is the frequency of the 
wave. The linear theory model for the interaction of par- 
ticles and waves is only valid for waves that are relatively 
long-lived; qualitatively, wave decoherence will cause the 
delta-function resonance to broaden, with more particles 
able to approximately satisfy the resonance condition for 
long enough to experience significant acceleration. 

Because TTD arises from the mirror force, it will be- 
come negligible as /i — > 0. In this limit (and in the high- 
limit, as we will show) the most important mecha- 
nism for changes in parallel velocity is what we r efer to 
as Fermi Type-B (FTB) acceleration (|Fermilll949[ ). Con- 
sider a particle spiralling along a magnetic field line with 
some curvature (see Fig. [I]). In the frame of the field 
line, the particle has constant energy, because magnetic 
fields do no work. In the frame of the bulk plasma (i.e., 
the frame in which the average fluid momentum is zero) , 
however, the particle may gain or lose energy. This can 
be seen straightforwardly with a Galilean transform (in 
the non-relativistic case) from the field line frame to the 
plasma frame. Qualitatively, FTB describes charged par- 
ticles as beads sliding along moving wires. These stochas- 
tic interactions will also cause diffusion in velocity space 
independent of [i. 

As noted above, we do not consider particles with 
Q ~ u). Particles with u) — k\\v\\ = nfl, where n is a 
non-zero integer, are "cyclotron-resonant" and can thus 
experience violation of /i-conservation and strong per- 
pendicular heating. However, due to the anisotropy of 
the st rong MHD turbulent cascade (iGoldreich fc S ridhar 
1995), MHD turbulence with a substantial inertial range 
may transfer its energy into a kinetic Alfven wave cas- 
cade on scales of k±, max ^ tT , where is the proton 
gyroradius. At this scale in the solar wind, the maxi- 
mum gyrofrequency that can be cyclotron-resonant for 
thermal particles is <C Cl p , the proton gyrofrequency (see 
ILehe et all [20091 ). As a result, any cyclotron heating of 
protons or electrons occurs (if it occurs at all) on scales 
below the MHD casc ade that we consider here. However , 
there is recent work ([Cranmer fc van Ballego oiienll2012h 
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suggesting that mode coupling between the isotropic fast 
mode cascade and Alfven modes could supply enough 
high-fcii power in the solar wind to heat protons through 
the cyclotron resonance. 

Our test particle calculations presented in Sj3] do not 
include diffusion due to parallel electric fields (Landau 
damping). The reason is that our simulations are per- 
formed in the ideal MHD limit. Physically, sound waves 
in a collisionless plasma do generate significant parallel 
electric fields - this effect is potentially important for 
the fast mode in fj > 1 MHD turbulence and for the slow 
mode in j3 < 1 turbulence, but is not captured in our test 
particle calculations. Additionally, the ideal MHD limit 
does not allow us to capture heating and acceleration 
by magnetic reconnection, which is likely an important 
mechanism in many as trophysical environment s includ- 
ing solar flares (see e.g. IDrake et aLl feOOG. 200l . 

3. ANALYTIC TRANSPORT PROPERTIES 

In this section, we provide an analytic derivation of 
the evolution of a distribution of test particles in veloc- 
ity space as they interact with Alfvenic MHD turbulence. 
We focus on the diffusive evolution in vu , which domi- 
nates over any changes in v±, as argued in <|2] 



3.1. Transport in v± 

The magnetic moment /j = Vj_ / B is approximately a 
conserved quantity in our test particle calculations de- 
scribed in J5] (for more discussion, see 

If /i is conserved, then any change in v± results from a 
change in the local value of the magnetic field. A distri- 
bution of particles that initially all have the same value of 
vj_ , randomly initialized in space throughout a turbulent 
plasma, will quickly broaden in v± over short time scales, 
and then reach a saturated width in v± as the particles 
statistically sample all of the fluctuations in B. Taking 
a differential of v± = ^/JlB, we find that final width of 
the distribution in v± will be approximately given by 
8v±/v± ~ Bl/2Bq, where Bj, is the rms deviation of B 
from the mean. We can understand this more simply by 
noting that if a particle is initialized where B is higher 
than J5o, then /x for that particle will be lower than the 
average. Thus the fractional width of the initial distri- 
bution in fi will be the same as the width of the initial 
distribution in B (in the Bl/Bq <C 1 limit). After some 
time passes and the test particles statistically sample the 
turbulence, a particle with a smaller fi is not preferen- 
tially likely to be in a location of larger magnetic field, 
so this particle will have a smaller time-averaged v±_ . To 
the extent that fj, is well-conserved, there will be little 
diffusion in v± after this initial re-adjustment. 

3.2. Transport in v\\ 

Diffusion in parallel velo city in our sim ulations comes 
from two sources, which iFermil (|1949f ) referred to as 
"Type A" and "Type B" interactions. Type A refers to 
acceleration by magnetic mirror forces. In Fermi's orig- 
inal conception, charged particles reflected off of mag- 
netic inhomogeneities (clouds) due to the fiVB force. In 
this work, the scattering centers are compressive MHD 
modes, so that Type A acceleration and TTD are effec- 
tively the same. Type B refers to the acceleration of 
a particle in the rest frame of the bulk plasma, due to 
following a curved, moving field line. 

3.2.1. TTD and Resonance Broadening 

A particle moving in compressible MHD turbulence 
will be randomly accelerated and decelerated by mirror 
forces, given by /zV B. This will lead to diffusion in par- 
allel velocity given by 

where /(x,v) is the velocity-space distribution function 
(we will assume a uniform spatial distribution of test 
pa rticles t h rough out). We adopt the fo rmal approach 
of iDupred (|1966l) and iWeinstockl (|1969T ) and write the 
corresponding diffusion coefficient as 

/>oo 

£|l=/! 2 / dt <V||B(r ,0)V||B(r(t),t)>, (2) 
Jo 

where () indicates an ensemble average, Vm means the 
gradient in the parallel (to the magnetic field) direction, 
and r(t) describes a particle's trajectory. Introducing 
a Fourier transform of the magnetic field, we can then 
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write Equation [5] as 

D " = \w J d3k k pM R( V> ^ 

where I(k) is the power spectrum of magnetic field fluc- 
tuations, and R(k) is a resonance function. If we consider 
only free-streaming trajectories (ignoring any variations 
in «||), as in linear theory, the resonance function is given 

by 

poo 

i?(k)=Re/ dt e'Kk)-^!!)*. (4) 

This resonance function becomes the standard S(v^ ±v p ) 
for interactions with waves satisfying the linear disper- 
sion relation uj(k) — ±\k\\\v p . The n^7\\B force re- 
quires compressive fluctuations, and thus will stem from 
slow or fast waves in MHD turbulence, rather than 
Alfven waves. However, the slow mode cascade pas- 
sively adopts the anisotropic G S power spectrum, wit h 
fej_ ^ fc|| in the inertial range (|Cho fc Lazarian|[200l . 
In this regime, the slow mode dispersion relation be- 
comes uj — k\\c s VAj V ' c 2 + v\ (|Chandranll2003t ). limiting 
to uj = k\\VA for $ S> 1 and uj — k\\c s for (3 <C 1. For 
analytical simplicity, we use a generic Alfvenic dispersion 
relation, uj = k^v p , where v p is the appropriate phase ve- 
locity. The specification of the results presented here to 
various slow mode regimes is straightforward. 

Deviation from free-streaming tr ajectories modifie s the 
resonance function in Equation @] (Weinstock 1969]); the 
resonance function becomes 

/>oo 

i?(k)=Re/ dt e »("(k)-« l |fc||)t+^[ l k.,5r(t)] 2 > ; (5) 
Jo 

where Sr(t) describes the deviation of a particle's trajac- 
tory, and we have assumed the mean deviation (Sr(t)) = 
0. Formally, this result requires kT B ^\\B <C 1, where 
tb is the correlation time of the stochastic magnetic field 
fluctuations on a given scale k, which will not hold for all 
scales or for all test particles we consider. Furthermore, 
the exact value of (\Sr(t)\ 2 ) depends on the diffusion coef- 
ficient that we are trying to calculate, and would require 
a recursive approach to calculating Dm and D±. Thus, 
for analytical simplicity, we use Equation [5] only to mo- 
tivate a phenomenological modification to the resonance 
function. 

To estimate an appropriate modification, consider that 
one source of deviation from free-streaming trajectories 
is the fact that high-/x particles are tied to magnetic 
field lines that can wander in the perpendicular direction. 
This process will be essentially stochastic and therefore 
diffusive. The typical step length will be given by vitb, 
where vi is the scale-dependent rms turbulent velocity. 
Thus: 

<<5ri(t)) ~ ( Vi tb) 2 (J^j ~ vfr B t. (6) 

Throughout the main body of the paper we fo- 
cus on the strong anisot ropic turbulence model of 
iGoldreich fc Sridharl (|1995l henceforth GS), though we 
con sider an extension to weak turbulence in Appendix 
IBl 5 ! In the GS model, the correlation time of fluctua- 
tions is given by Tg 1 ~ ui n \ — (kj_L) 2 / 3 VL/L, where vl is 

5 Note that since we drive our test particle simulations in J6] 



the rms turbulent velocity on the driving scale L. This 
timescale may be thought of as the "lifetime" of waves in 
turbulence. Thus our estimate for the field line random- 
walk deviation becomes 

([k • 5r(t)} 2 ) ~ v L (k ± L) 2 / 3 t/L = cj nl t. (7) 

The equality above suggests that we could have reached 
the same conclusion via the somewhat different approach 
of directly modifying the time dependence of the wave 
modes in our turbulence: e lu)t — > e 4W * _ " nl * (defined for 
t > 0). This alternative approach directly describes the 
modes themselves as decohering on a timescale uj~\ . We 
will refer to this as the exponential decorrelation model. 

To account for the uncertainty in these estimates, we 
replace the factor of 1/2 in Equation [5] with 7, a di- 
mensionless order-unity constant. In this section we also 
choose to instead use ui n \ = {2nv A /L)(k ± L/2n) 2 / 3 , the 
generic "strong turbulence" expression for the non-linear 
turnover timejj This is for the sake of generality: for 
MHD turbulence with a significant inertial range, even if 
the turbulence is weak on the outer scales, the nonlinear- 
ity of the turbulence increases on smaller scales, eventu- 
ally approaching critical balance. When critical balance 
is reached, the remainder of the cascade will be in the 
strong regime. Furthermore, the turbulence is weak on 
the outer scales, where the turbulence is being driven. 
Thus, the details of the driving are more likely to be im- 
portant for the turbulence statistics in the weak regime. 
We present analytic calculations of resonance broadening 
in weak turbulence in Appendix [Bl 

Given Equation^ we can perform the integral in Equa- 
tion [5] to calculate the exponential resonance function, 

d n , 7 ^(fc ± £) 2 /3 (27r) i/3 

"-PW 7 2 u 2 i (2 7r )2/3( fc±i )4/3/ jL + fc 2 L(u|| ± Vp) 2 > 

where we have used the dispersion relation, ui(k) = 
±Uj,|fcii|. We see that the resonance function is still 
peaked at |um| = v p , but becomes a Lorentzian in U11 
rather than a delta-function. 

To calculate the parallel diffusion coefficient Dm asso- 
ciated with this broadened resonance, we perform the 
integral in Equation [J2 using a power spectrum of the 
GS form, 

J(k) oc vlk^g (^Pj , (9) 

normalized such that v\/2 = J d 3 k I (k) / (2tt) 3 . In this 
expression, we will treat g(x) as a step function, equal 
to 1 if \x\ < 1, and otherwise, which accounts for the 
fact that power only resides in ku -C k± in the iner- 

at sub-Alfvenic velocities on the outer scale, we are in fact in 
the weak/intcrmcdiatc turbulence regime on size scales of l± > 
l c = LM'j^, where L is the outer (driving) scale of the turbu- 
lence and M\ refers to the Alfvenic M ach number at L (see e.g. 
IGoldreich fc Sridharl[T997l: ICJaltier et alj|200a) . We discuss an ex- 
tension of the model in this section which incorporates weak tur- 
bulence in Appendix [51 

6 If one considers a simulation with fixed k and va (i.e., a fixed 
simulation box), decreasing e (and therefore v^) corresponds to 
weaker turbulence, since w n i < to ~ k^vj^. 
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tial range of the cascadeQ The turbulence model de- 
scribe d here is broadly supported by numerical si mula- 
tions (iMaron fc Goldreichf 120011: iCho fc Lazarianl l2002t 
iBeresnvald I2012D. H oweve r, it is not univers a lly ac - 
cept ed; see IBoldvrevI (l200ll. IPerez fc Boldvrevl (|2008l) 



and IGrappin fc Mullerl (|201oj) for discussions of possible 
shortcomings of this model. Substituting Equation[5]into 
Equation [3] and performing the integral leads to 



£>H = dGv] 



1 ( 1 

1 arctanu_ 



1 



■til 



1 



— =- 1 arctanwj 



(10) 



where u± = (uii ± v p )/ r yvA, C\ is a dimensionless con- 
stant, G = v\ 7rlog (L / L m i n ) / (6 j L v A ) is a function 
which absorbs normalization constants, and -L m j n is the 
smallest resolvable length scale on the grid. We leave the 
normalization constant C\ unspecified at the moment; if 
our previous calculations were exact, then C\ would be 
equal to 1. We will calibrate this value against our test 
particle simulations in SjHl Equation [TU] has the limiting 
values 

const if |f|| | <C v p 



D 



■x 



v pI v \\ 



I if|»||l>« 



v 



(11) 



Again, we see that the delta-function resonance predicted 
by linear theory is substantially broadened, so that all 



particles with 



< 



va couple equally well to the tur- 



bulence; in addition, high-velocity particles can also in- 
teract with the turbulence via the TTD resonance. 

It is not entirely clear on theoretical grounds that the 
exponential decorrelation of the preceding discussion is 
the correct or the only model for resonance broadening. 
Thus we consider also a simple alternative, which we re- 
fer to as a Gaussian decorrelation model. We replace 
e -7w nJ t m Equation IH1 with e - ^ 7 ^" 1 ^ , where 7 has quali- 
tatively the same physical interpretation as in the expo- 
nential case. We argue that this functional form for the 
wave decoherence is more physically motivated as it has 
smooth derivatives as t — > 0. In this case, we find the 
Gaussian resonance function, 



Ra 



exp 



jVA(4:k±L) 2 / 3 

k\L 2 {v\\±v p ) 2 



4 7 2 V 2(27r) 2 / 3 (fcj_L) 4 / 3 



(12) 



so that the (5-function becomes a Gaussian resonance. 
Again we use Equation [3] to find 



D« = C 2 Gvi 



+ 



1 



7r erf 



U- 



7r erf 



— u-v^exp 



(13) 



7 The step function approximation for g(x) is for analytic con- 
venience. Physically, the cutoff in power in the fcy direction is 
unlikely to be quite so sharp. We have confirmed using numeri- 
cal calculations that a cutoff in fcii that is, e.g., exponential, rather 
than a step function, produces quantitatively similar results for the 
diffusion coefficients of interest. 
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Fig. 2. — Analytic models of the effects of resonance broaden- 
ing on the parallel velocity diffusion produced by the interaction 
between particles and strong MHD turbulence. We show both ex- 
ponential (black lines) and Gaussian (blue) models of resonance 
broadening (see i]3.2.1l L with different values of 7 indicated by 
linestyles. Both are arbitrarily normalized, with the exponen- 
tial models artificially normalized to a higher value for clarity. 
v p = 0.81 c s is chosen, consistent with /3 = 1 slow modes. The 
dimensionless parameter gamma controls the magnitude of the 
broadening with gamma <C 1 approaching the linear theory pre- 
diction of a delta function at Vn = v p . The exponential resonance 

function gives ZJy <x i>~ 2 in the high-Dy regime, while the Gaussian 
model gives D« <x v7 3 . 



where erf(x) is the error function and C2 is again a 
dimensionless normalization constant to be calibrated. 
This has the limiting values 



const if |u|| I <C v p 

X " " 3 /l-l 3 if |v||| 



(14) 



In Figure [2] we plot several representative examples 
of Du with arbitrary normalization. The dimensionless 
parameter 7 controls the "peakiness" of the resonance. 
Note that in all cases D\\ declines steeply above v p , which 
implies that TTD heating of very fast particles is ineffi- 
cient. TTD acts primarily on particles in the bulk of the 
plasma, near the (linear) resonance. 

For completeness, we also calculate the diffusion coeffi- 
cient resulting from the linear theory delta-function reso- 
nance, with u} n \ = 0. Is this case, the resonance function 
becomes 



R(k) 



\k 



jS(v\\ ± v p ) 



(15) 



Once again, we apply the GS power spectrum to the 
resonant diffusion coefficient of Equation [3] to find 



D 



24L 



va 



v ± log 



S(v l{ ±v p ). (16) 



3.2.2. Resonance broadening of other modes 

Our discussion up to this point has focused on slow 
modes with j3 > 1 for the sake of analytical simplicity. 
These results are, however, relatively easy to general- 
ize, and can be applied to other wave modes and plasma 
parameter regimes. For example, analytically account- 
ing for fast modes is not difficult. In Figure [3l we plot 
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Fig. 3. — Fast mode contribution to resonance-broadened TTD 
assuming either an exponential or Gaussian for m for the decorre- 
lation (arbitrary normalization, with 7 = 1; see i]3.2,2l for more de- 
tails on the calculation). For comparison, we include also fast mode 
TTD with a delta-function resonance. The weak non-linearity of 
fast modes in MHD turbulence leads to a sharp resonance, very 
close to the linear theory result. 



numerically-calculated resonance-broadened TTD coeffi- 
cients associated with fast modes in /3 < 1 turbulence. 
We include for comparison the diffusion coefficient pre- 
dicted by a linear theory delta function resonance. (The 
linear theory diffusion coefficient in this case is not itself a 
delta-function due to the fast mode resonance condition, 
kvA — Several modifications to the derivation in 

13.2.11 are required for the resonance-broadening calcula- 
tion. We use the fast mode dispersion relation, uj — t^fc, 
and we assume an isotropic fast mode power spectrum, 
I(k) oc fc~ 7 / 2 (e.g. iCho fc Lazarianll200l . More im- 
portantly, the non-linear decorrelation frequency for fast 
modes, vik (where vi is the turbulence velocity on scale I), 
is much smaller than the corresponding linear frequency 
va fc, implying that fast modes decorrelate much more 
slowly (in turbulence) than Alfven and slow waves. Fur- 
thermore, this non-linearity becomes weaker on smaller 
scales. Thus, TTD with fast modes will be much less 
broadened. Figure [3] shows that this results in a veloc- 
ity diffusion coefficient much more peaked near um ~ va, 
closer to the linear theory result. The broadened reso- 
nance does still result in a power-law diffusion coefficient 
at high vh , with the same power-law indices as in the slow 
mode case considered above (-2 and -3 for the exponential 
and Gaussian decorrelation models, respectively). How- 
ever, the distinction between the exponential and Gaus- 
sian models only becomes apparent at very high vu . This 
high-uii tail proves to be the most important feature of 
the broadened resonance for calculating heating rates at 
high v s . 

Similar modifications to the calculation in £|3.2.1l would 
need to be made in other regimes. For instance, for slow 
modes with /3 < 1, Equation|H]would need to be modified 
by a multiplicative factor ~ (3/(1 + /3) to account for the 
decreasing magnetic compression of slow modes in this 
regime. Additionally, one would instead use a dispersion 
relation uj — ±c s |fcii|. However, we anticipate that the 
general functional form of the diffusion coefficient is sim- 
ilar to those shown in Figures [2] and [3] in these different 



regimes. We will use this generality of the resonance- 
broadened diffusion coefficients to interpret our test par- 
ticle simulation results in $6l 

3.2.3. Type B diffusion: the n — ¥ limit 

In the fi — > limit, magnetic mirror forces become neg- 
ligible and all diffusion in parallel velocity is due to Fermi 
Type B interactions, resulting from the tying of particles 
to moving, curved field lines. In our non-relativistic case, 
this change in parallel velocity in one coherent interac- 
tion with a curved field line will typically be of the order 
5vu ~ VLSmO ~ vl(SB / Bq), where 9 is the opening an- 
gle of the magnetic field line, as illustrated in Figure [TJ 
These interactions will be stochastic, and we estimate a 
parallel velocity diffusion coefficient by D\\ ~ Sv 2 /t COIV , 
where t corr is the typical time over which a particle expe- 
riences correlated field line motion. For a particle with 
\ v \\ \ *C vl,va, the decoherence time of interactions will 
be determined by the outer-scale fluid motions. The el- 
ements of field line curvature which provide FTB diffu- 
sion may be thought of as essentially Alfvenic fluctua- 
tions, because they are most effective when <5B _L Bo- 
For strong Alfvenic turbulence, the decoherence time of 
a wave-particle interaction at the outer scale will be pro- 
vided by a combination of two effects: linear propagation 
and non- linear distortion (eddy turnover). For weaker 
turbulence, the linear propagation of fluctuations will 
control the decorrelation of wave-particle interactions. 
Furthermore, at any instant in time, the outer-scale fluc- 
tuations have correlation lengths ~L, because the turbu- 
lence is driven on this scale. Thus in either case (strong 
or weak), a good estimate of the wave-particle correla- 
tion time is t corr ~ L/va, the outer-scale wave crossing 
time. This allows us to estimate the diffusion coefficient 

'SB"" 2 
Bo. 



Dn 



L/va 



Lv A 



(17) 



where the second equality follows from assuming that the 
typical magnetic field perturbation at the outer scale is 
of order 5B / Bq ~ vl/va, which will be true for Alfvenic 
turbulence. We could choose to express this in terms 
of e, the turbulence cascade rate, which in Kolmogorov- 
or GS-like turbulence scales as e ~ v\/L, However, the 
large-scale eddies are those most effective at FTB accel- 
eration. The strong turbulence scalings are least likely 
to be applicable on these large scales, and so we leave 
the FTB diffusion coefficients explicitly in terms of vl- 

On the other hand, particles with \v\\ \ 3> va are essen- 
tially interacting with a static snapshot of turbulence, 
and so i corr ~ L/|ui||, the particle crossing time of the 
outer-scale correlation length. This implies 
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5BY 
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(18) 



We choose a functional form for Dn that asymptotes 
to the scalings in Equations [T7] and [18] in the low- and 
high- velocity regimes and varies smoothly between these 
limits: 



FTB 



= c. 



Lv A 



C 4 



(19) 



C3 and C4 are dimensionless parameters which we will 
calibrate against test particle simulations in <J6] 
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3.2.4. Phase-decorrelation broadening 

The discussion in §3.2.11 focuses on the phenomeno- 
logical idea of wave-particle phase decorrelation as a re- 
sult of the decay of the resonant wave. However, the 
changes in parallel velocity which the test particle ex- 
periences in the wave-particle interaction will also lead 
to phase decorrelation. As a particle's parallel velocity 
changes, its position z along the magnetic field changes 
relative to a ballistic trajectory with z\, = zq + v\\ t. 
The difference Sz — z — Zb is given by Sz — J Svudt, 
where 5vu is the change in particle velocity resulting from 
the acceleration process. When Sz ~ 1/ku, the particle 
has moved completely out of phase with the wave. For 
parallel velocities which change diffusively, we may esti- 
mate that t he rm s change in parallel velocity is given by 
fiv\\,ims ~ v/^p- Estimating Sz ~ Sv\\ trms t, the typical 
extra random phase between wave and particle will be 
4> ~ k\\D^ 2 t i l 2 oj (i/i ph ) 3 / 2 , where the final approxi- 
mate equality is simply a definition of t p h- One could 
now use this additional decorrelation in an extension of 
the broadened resonance of Equation [5] However, this 
introduces a recursive dependence of Dm on itself, seem- 
ing to limit the analytical tractability of this approach. 
For simplicity, we do not include this effect in our an- 
alytical model. However, we do consider this effect in 
interpreting our test particle results. 

Using instead <fi 2 ~ k 2 D^t 3 in the resonance broad- 
ening calcu l ation is formally similar t o Equ ation 61 of 
IWeinstocM (|1969Q . lYan fc Lazarianl (|2008Q also used 
a similar approach in modelling resonance broadening. 
The equivalent of i p h which they calculate may be found 
by using a ballistic approximation for the particle devia- 
tion, rather than a diffusive one, so that 5v\\ rms ~ /iV Bt 
(though this deviation is still treated as effectively ran- 
dom over a distribution of particles, in that it may be 
parallel or anti-parallel to the mean magnetic field), in 

1 1/2 

which case the effective t ph ~ k\\v±M A ' . This ballis- 
tic assumption may be more be appropriate for parti- 
cle transport at early times. (This scaling assumes that 
5|-Bj cx Ma, which is not true for pure Alfven waves but 
is the case in MHD turbulence with a significant com- 
ponent of compressive energy.) This i p h is essentially 
identical to the bounce time for a particle of magnetic 
moment /i in a magnetic wave of amplitude MaBq and 
wavelength ~ I/&11. 

4. HEATING OF A THERMAL DISTRIBUTION 

If the evolution of a distribution of test particles satis- 
fies a diffusion equation (as in eqn. Q]), we may multiply 
both sides of this evolution equation by m s v 2 /2 and in- 
tegrate over all velocities to find the volumetric heating 
rate of the particles, 
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(20) 



where T s is the temperature of the particle species under 
consideration, fee is Boltzmann's constant, and the dis- 
tribution is normalized such that J d 3 vf s (v) — n 8 , the 
spatial density of particles. For simplicity, we assume 
that the distributions are Maxwellian, with a thermal 



velocity v s — -y/fcsT s /m s . We will treat T s — T as a 
constant when we compare heating of different particle 
species, appropriate for species in temperature equilib- 
rium. 

4.1. FTB Heating 

The parameterization of FTB diffusion in Equation[T9l 
as well as the assumption of a thermal distribution, al- 
lows us to calculate the FTB heating rate from Equation 



E s ,ftb — ntcstksT 
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(21) 



n tcs t refers to the number density of the test particles. 
The v~ 2 dependence of this expression at small v s will 
cause FTB to dominate the heating of low thermal ve- 
locity particles (v s <C c s , as for e.g. minor ions). 

4.2. TTD Heating: Linear Theory 

The linear theory (LT) diffusion coefficient for slow 
modes is given by 

Dm oc i^l^C^H w p)' 
which gives a heating rate 



(22) 



E sX T oc Tv s 1 exp I --^2 



(23) 



For v s ^> v p and species at roughly the same tempera- 
ture, this heating rate scales as v^ 1 . In the next section, 
we find that resonance broadening in general implies a 
shallower dependence on v s for the asymptotic heating 
rate of high-w s particles. 

4.3. TTD Heating: Resonance Broadening 

Substituting the exponential and gaussian resonance 
broadening expressions into the heating rate integral in 
Equation [201 does not lead to a simple analytic integral, 
and so in our comparison to our test particle simulations 
we will numerically evaluate Equation 1201 However, as 
we will see, TTD is the dominant contribution to the 
heating for high velocity particles. Thus we may gain 
some insight by considering the heating in the v\\ S> v p 
limit. 

In particular, the exponential decorrelation function 
gives D|| cx w l w ||^ 2 for «|| ^> v p . This implies a high-w s 
heating rate of 

E s ^ttd cx T, (24) 

independent of v s . Thus, the exponential model leads 
to more effective heating for high-velocity particles (e.g. 
electrons) than linear theory. The Gaussian decorrela- 
tion function gives Dm oc v^|i>|j|~ 3 for U11 3> v p , on the 
other hand, which implies a high-t; s heating rate of 



E s ,ttd cx 



T\n{v s ) 



(25) 



which has a scaling intermediate between the LT and 
exponential cases, although closer to the linear theory 
result given that the only difference is the weakly varying 
\a{v s ) factor. 



5. NUMERICAL METHODS 

Our simulations consist of collisionless test parti- 
cles evolving in isothermal, subsonic MHD turbulence. 
Our compu t ation al approach is quite similar to that of 
iLehe et all ()2009[ ). apart from two important changes 
noted below. We present a summary of our methods 
here; more detail may be found in the earlier paper. 

5.1. The MHD Integrator 

We use the Athena MHD code of lStone et all (|2008fl to 
evolve the turbulence on a 3D Cartesian grid with peri- 
odic boundary conditions. The grid is initialized with a 
uniform background magnetic field Bq in the x-direction, 
with the velocity set to zero everywhere. The initial mag- 
nitude of B is set by our choice of /3 = pc 2 s / (B 2 / '8ir) , 
where p is the fluid density and c s is the sound speed. 

We then inject kinetic energy by providing "kicks" 
to the velocity field, in a method similar to that of 
iLemaster fc Stonel (|2009D . At each timestep, we generate 
a velocity perturbation <Jv(k) with random amplitudes 
in Fourier space in the range of2x^<fc<4x^, 
normalized by a decreasing power law in k, so that the 
majority of the driving power enters on the largest scale 
of L/2. We also remove modes with |fcii| < 2 x to 
avoid parallel correlation lengths longer than ~ L/2. We 
enforce <5v(k) • k = 0, so that our velocity field is di- 
vergenceless, to minimize the excitation of compressible 
modes (see the discussion at the end of this subsection 
for more detail on the decomposition of the turbulence 
in MHD modes). We then normalize <5v(k) so that the 
net energy input into the turbulence is given by e. 

We ensure that the kicks are time-correlated by imple- 
ment ing an Ornstein-Uhlenbeck (OU) process (|Bartoschl 
[200lh . given by 

5v(k, t + dt) = /<Sv(k, t) + ^/l-f 2 Sv'(k), (26) 

which has an autocorrelation time (assuming the contin- 
uous limit, dt — > 0) given by 

(*v(Mi) -5v(k,h)) = {[5v k {t)] 2 )e^-^l^" (27) 

where / = e~ d */* corl , t COII is the correlation time of the 
driving, dt is the timestep of the driving routine, and <5v' 
is a new random field generated by the process in the 
previous paragraph. We choose to drive on every MHD 
timestep. The OU process is simply a mean-reverting 
random-walk. Note that in order for Equation [27] to 
properly describe the driving statistics, the initial kick 
<5v(k, t — 0) must be drawn from the same random dis- 
tribution as the subsequent <5v (k). 

Time-correlated driving is critical for two reasons. 
First, any process that drives turbulence on large scales 
will be correlated on some typical timescale depending on 
the underlying physics of the driving process, rather than 
pure white noise. Thus a time-correlated driving scheme 
is more representative of the underlying physics of the 
turbulence. More pragmatically, evolving our test par- 
ticles in turbulence with <5(i)-correlated driving leads to 
unphysical acceleration of high gyrofrequency particles, 
because of the high frequenc y power pre sent in the tur- 
bulent driving. To avoid this. ILehe et al.l restricted their 
analysis to test particles interacting with decaying (non- 
driven) turbulence. Driving via the OU process allows 



us to consider particles evolving in saturated turbulence 
over arbitrary lengths of time. 

On physical grounds, we choose to apply a correlation 
time of order L/vl, the eddy turnover time on the outer 
scale of the turbulence; see ^6.61 for a fuller investigation 
of the dependence of particle heating on the correlation 
time. Additionally, we must choose < corr ^> u^, where 
w max is the maximum wave mode frequency resolvable 
in the MHD simulations. Smaller values of t corr imply 
essentially uncorrelated driving and lead to unphysical 
heating through a resonance with the MHD timestep. 

For a simulation with periodic boundary conditions 
and a velocity field driven on the size of the domain, a 
particle with arbitrarily high velocity effectively encoun- 
ters the same eddy repeatedly, as it crosses the box many 
times before the eddy decorrelates. This is unphysical, 
and thus we choose a fiducial volume for our simulations 
of {16L, 2L, 2L}, so that the box is elongated in the di- 
rection parallel to Bo, and there are approximately 32 
uncorrelated eddies along the length of the box. If wc 
use instead a cubical box of side length 2L, we find D|| is 
unphysically affected by box-crossing for particles with 
velocities vn > 10 c s . Extending the box to a size of 
length 16L in the parallel direction allows us accurately 
evolve particles with velocities up to ^50c s . This is 
particularly important for studying the evolution of elec- 
trons, corresponding to our high- velocity particles. 

The parallel extension comes at the cost of decreased 
resolution at the smallest scales. However, FTB acceler- 
ation is dominated by the largest eddies, and therefore 
accurately capturing smaller eddies is irrelevant to zeroth 
order. Similarly, slow-mode TTD has only a logarithmic 
dependence on the length of the inertial range. Thus we 
choose to focus computational resources on the larger- 
scale eddies. 

We choose p = c s = L = 1, but the results of our 
simulations can be applied to different physical systems 
by scaling them with appropriate combinations of p, c s , 
and L. Thus our turbulence is controlled by three pa- 
rameters: the specific energy input rate e, in units of 
c s /L; the ratio of plasma to magnetic pressure (3; and 
the correlation time t CO rrj in units of L/c s . 

We note for reference that we have applied the approx - 
imate, Fourier-space method of lCho fc Lazarianl (2003) 
to decompose the turbulent kinetic energy in our simu- 
lations into Alfvenic, slow, and fast mode components. 
Across the range of driving rates in our simulations at 
fixed /3 = 1, roughly 45% of the kinetic energy is in 
Alfvenic modes and 45% is in slow modes. At lower /3, 
an increasing fraction of the total energy is in slow modes, 
up to 60% for f3 = 0.1, while Alfvenic modes lose a corre- 
sponding fraction. Fast modes never comprise more than 
^5% percent of the kinetic energy, and a similarly small 
fraction belongs to motions with — which cannot be 
identified with any MHD wave mode, corresponding to 
interchange modes. Changing the correlation time also 
results in somewhat different 2D power spectra; specifi- 
cally, longer correlation times appear to frequency- match 
onto low-fcy modes, so that when the turbulence saturates 
there is power in modes which are not directly driven. We 
discuss this further in Appendix [Cl 

5.2. Particle Integration 
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Our par t icle i ntegration methods are described in 
iLehe et all (|2009ft . Once the turbulence reaches a fully- 
saturated state, particles are evolved according to the 
Lorentz force. We describe particles by their charge-to- 
mass ratio, expressed in the form of the mean gyrofre- 
quency Qq = qB^/mc. The actual gyrofrequency of a 
particle will vary according to the local value of B, but 
in subsonic turbulence, SB /Bo < 1, so variations in Q 
are not large. Our simulations use ideal MHD, with the 
resistivity r\ set to zero, so the turbulent dissipation is nu- 
merical. Thus the electric field is given byE = — vxB/c, 
where v is the fluid velocity. 

We integrate the particles with the iBorisI (|197O0 im- 
plicit particle pusher. This method is symplectic and 
symmetric in time, and conserves energy and adiabatic 
invariants to machine precision in simulations with con- 
stant fields in space and time. We choose a timestep 
much smaller than the gyroperiod of the particle. We 
interpolate the MHD fields on the grid to their value 
at the particle's location using th e Triangular-Shaped 
Cloud (jHocknev fc Eastwoodlll98ll ) method in space and 
time, while ensuring that the interpolation does not in- 
troduce spurious parallel electric fields (E\\ = 0). 

We initialize the particles randomly over the simula- 
tion volume, assigning them a v± and Vu , where these 
are measured perpendicular and parallel to the local mag- 
netic field0 The perpendicular motion of a particle is the 
superposition of the fast gyration around B and a slowly- 
varying drift velocity, v±_ = |v to t - Yd|, where v to t is 
the total perpendicular velocity of the particle and V£> is 
the drift velocity. We thus require knowledge of the local 
drift velocity to accurately assign v± . We account for the 
E x B drift, VB drift, curvature drift, and the polariza- 
tion drift. These latter three drifts are typically smaller 
by a factor of uj/Q,, where w _1 is the timescale of a fluc- 
tuation, so one might naively expect them to be small in 
our simulations. However, the curvature drift, approxi- 
mately given as vr> ~ v 2 /QR c (jHazeltine fc Waelbroeckl 

1998), where R c is the local radius of curvature of the 
magnetic field line, can become important at high v\\ . 

For many of our simulations, we initialize a distribution 
of particles with fixed v± and a logarithmically-binned 
distribution in v\\ (or vice- versa), to isolate the effects of 
one variable. In other cases, we initialize particles ac- 
cording to an isotropic Maxwell-Boltzmann distribution: 



TABLE 1 

Summary of fiducial simulation parameters 



exp 



II 



2v 2 



(28) 



where v s is the typical thermal velocity of the distribu- 
tion. We summarize fiducial parameters for our simu- 
lations in Table [TJ and explicitly note elsewhere when 
different parameters are used. 

6. TEST PARTICLE DIFFUSION IN SIMULATIONS 

We initialize a distribution of particles in vi\, v± or 
v, the magnitude of the velocity. For particles initially 
within a given bin in v±_ and vu , we calculate diffusion 
coefficients according to the formal definition 

St— >oo ZOt 

8 Our TSC interpolation scheme means that the local magnetic 
field is measured on approximately the grid scale. 



Parameter 



Value 



P 

Cs 

L 

Resolution 
Volume 
i (c 3 e /L) 


-^particles 

n (c s /L) 

icorr (_L / Cs) 



i 
i 
l 

1024 X 128 2 
16L x (2L) 2 
0.1^ 
1~ 

2 n x 10 3 ~ 2 x 10 5 
10 5 
1.5 



a This produces a sonic Mach number of ~ 0.35. 

where A is the diffusing quantity. We are typically in- 
terested in calculating Du , the parallel velocity diffusion 
coefficient. Our fiducial set of results are for /3 = 1 tur- 
bulence with e = 0.1c^/L (this holds for Figures [4][9]) . 
We measure the diffusion coefficients over a time dura- 
tion from test particle initialization until after the initial 
ballistic behavior has become diffusive. For simulations 
presented here, this is typically between 0.1 and 0.75 
L/c s (with shorter durations for higher turbulent ampli- 
tudes). 

6.1. Non-conservation of fi? 

We assumed throughout our analytic calculation in Sj3] 
that \x is conserved. We do observe diffusive changes in /x 
throughout our simulations; i.e., /i is not in fact strictly 
conserved. However, the changes in fi we find do not sig- 
nificantly affect our parallel diffusion results or our in- 
terpretation of these results. In Figure IH we plot D^/fi 2 
for our fiducial simulation. All parallel velocities expe- 
rience some diffusive change in fi. However, this change 
is fractionally small until the highest vu . Furthermore, 
for particles with vu 3> v±, velocity diffusion is primar- 
ily due to FTB, which is independent of /x, so that /x 
changes significantly only in regimes where it is irrele- 
vant to the dynamics. Moreover, the diffusion coefficient 
for /i is fractionally much smaller than the corresponding 
diffusion coefficient for vu . Thus we are justified in using 
the approximation that /j, is conserved. 

6.2. Diffusion in vu 

Figure [5] shows our calculated D\\ for particles with 
v± = 0.3 c s (solid line). Our fiducial value of e is suf- 
ficiently small that [J.VB forces are almost negligible 
for small values of v± and therefore diffusion is dom- 
inated by FTB. (A similar run with v± — 0.1 c s , not 
shown, is essentially identical.) At low v\u the dif- 
fusion coefficient saturates to a constant value of or- 
der v\jLvA ~ 0.005 c^/L, consistent with Equation [T"7l 
(though a factor of ~ 4 smaller). At high v», the diffu- 
sion coefficient is proportional to v\\ , consistent with our 
analytic derivation in Equation 1181 

Figure [5] also shows that for particles with larger v± 
(larger /i), Du is significantly larger for particles with 
v \\ ^ fsiow, the phase velocity of slow modes. This is 
due to TTD, which increases in importance for larger 
fi. In particular, for larger v±, the TTD contribution 
manifests itself as an approximately constant Du for < 
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Fig. 4. — Diffusion coefficient for the magnetic moment fi, for 
several values of the initial perpendicular velocity. Diffusion in 
/i is fractionally much smaller than diffusion in tlii; furthermore, 
changes in »x are dominated by local variations in the magnetic 
field rather than the observed slow diffusion in fi. Thus in gen- 
eral we are justified in approximating the magnetic moment as a 
conserved quantity. 



^siow, and then as a smooth decrease in Dii for higher 
parallel velocities. This is consistent with the resonance 
broadening TTD models in M3.2.1I 

For particles with vu 3> fj_ beyond the linear resonance 
at Vu = Wsiow, FTB begins to again dominate the parallel 
velocity diffusion. The high vh scaling of Dn is oc vn, 
consistent with Equation [THJ The importance of FTB 
can also be seen by the fact that all of the curves in Fig 
5 are the same at high i>||, independent of v±. This is 
because FTB rather than TTD provides the dominant 
source of velocity diffusion at high vn, and the value of 
v± (//) is irrelevant to the efficiency of FTB diffusion. 

The broadened resonance in Figure [5] appears to move 
to the right for increasing v±. We argue that this is 
the result of phase-decorrelation broadening, discussed 
in £)3.2.4I We initialize delta-functions in vn, but as a 
result of the changes in vu caused by finite-amplitude 
turbulence, these bins quickly begin to sp read out. We 
may set the decorrelation time of £13.2.41 equal to lu 

fc||V| 



A)|| (v p — vh) to find the resulting broadening width 
Aw|| = v p — v\\ . The phase-decorrelation model in which 
the initial particle transport is ballistic in vu (as indeed 

1/2 

we observe at early times) predicts that Av\\ ~ v±, 
which is consistent with the test particle results in Figure 

El 

Figure [5] shows Dn as a function of v± for a distri- 
bution of particles with vtij = 0.2 c s . For high vj_, the 

diffusion coefficient scales like Vj_, consistent with the fi 2 
scaling in Equation 1101 This scaling may be understood 
by noting that D\\ oc a 2 oc (/iV£>) 2 , where a\\ is the in- 
stantaneous acceleration felt by a charged particle with 
magnetic moment /i. For smaller v±, as /i — > 0, the diffu- 
sion reaches the floor provided by the FTB mechanism. 
However, we note that at later times, we see shallower 
power laws in v±. We believe that this is also caused by 
phase-decorrelation effects due to finite changes in parti- 
cle uii discussed ^3.2.41 
Figure [7] presents a more quantitative comparison be- 



FlG. 5. — Parallel diffusion coefficients as a function of parallel 
velocity vn for several values of v± > i (in units of c s ) for our fiducial 
turbulence properties in Table [T] For the lowest ; = 0.3 c s 
(solid curve), diffusion is provided essentially entirely by the FTB 
mechanism, which provides an effective floor for Dy at low v±. 
For higher v±, resonance- broadened transit-time damping causes 
further diffusion, characterized by a constant Dy below the linear 
theory resonance at vn = v s \ ow ~ c s , and a smooth falloff above. 
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Fig. 6. — Parallel diffusion coefficient measured at early times as 
a function of perpendicular velocity for a distribution of particles 
with Dm = 0.2 c s (other properties of the simulation are summarized 

in Table [TJ. At high v±, Dy <x ujj this is due to transit-time 
damping. At low v± , the diffusion reaches a floor provided by the 
Fermi Type-B mechanism. At later times, we see shallower power 
laws in v± . We b elieve that this is caused by phase-decorrelation 
effects; see i|l2. 11 



tween our numerical test particle results and the analytic 
results for Dn derived in H3.2.1I The curve labeled by 
TTD refers to diffusion with a functional form provided 
by the Gaussian decorrelation prescription of Equation 
[T3l FTB refers to the sum of the contributions of eqs. [T71 
and llSI The normalization of the analytic diffusion coef- 
ficients are chosen by eye so as to provide the best match 
to the numerically determined diffusion coefficients. Our 
analytic model captures the qualitative character of the 
test particle results. 
Figure [8] shows the time dependence of our test particle 
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Fig. 7. — Comparison of numerically calculated diffusion coef- 
ficients to analytical scalings (using the Gaussian decorrelation 
model with 7 = 0.5), for a simulation with vj_ i = 3.0c s . The 
turbulence properties are given in Table [TJ TTD dominates at low 
dm, while FTB dominates at high vu. 
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Fig. 8. — Parallel velocity diffusion coefficient vs. v\\ , measured 
at several different times (in units of L/c s ) , for a simulation with 
e = O.OlCg/L, v± i = c s , and t CO rr = 3L/c s , approximately equal 
to the eddy turnover time. 

diffusion coefficients measured over different time base- 
lines, for a simulation with e = 0.01 c^/L and v± = c s . 
Perfectly overlying curves measured at different times 
would indicate perfectly diffusive behavior. We observe 
time dependence which is generally not perfectly diffu- 
sive, and is typically somewhat subdiffusive, in that e.g. 
((8v\\) 2 ) scales somewhat less than linearly with t. 

6.3. Approach to quasilinear theory? 

For smaller turbulence amplitudes, one might in prin- 
ciple expect the velocity diffusion coefficient to approach 
the sharp resonance of quasilinear theory, because the 
turbulence becomes increasingly weak on the outer scale. 
For example, in the simulation plotted in Figure [5J the 
amplitude of the turbulence is such that the turbulence 
is weak on the outer scale, with vl — 0.16 c s <C va- 
Our model of a weak turbulence cascade in Appendix [B] 



predicts a rather sharp resonance, for these parameters. 
However, no obvious resonance is present in the test par- 
ticle results shown in Figure [8J 

We believe that this is the result of phase-decorrelation 
broadening, as discussed in £13.2.41 This decorrelation ef- 
fect acts in addition to the primary wave-decay decor- 
relation model d iscuss ed in this paper. In turbulence, 
as discussed in ^3.2.41 there is a typical wave-particle 
phase decorrelation time given by t~^ ~ k 2 ^ 3 D^ 3 (if the 

particle velocity change is diffusive) or t~£ ~ k^v±M^/ 2 

(if ballistic), and we can use this to predict a broaden- 
ing width At)|| by equating the linear frequency with the 
decorrelation frequency, as in M6.21 

For the diffusive phase-decorrelation broadening, 
A«|| ~ (Dn/fcn) 1 / 3 . For the results in Figure |U this 
approach predicts Av» ~ 0.09 c s , evaluated at kit — 
4tt/L. This is well less than the measured broadening. 
On the other hand, the ballistic (bounce-time) phase- 

1 /2 

decorrelation broadening gives Av\\ ~ vxM% ~ 0.34 c s , 
which is consistent with the measured broadening to 
within a factor of ~ 2. We also note that the double- 
peaked features in Du in Figure [5] are similar to those 
plotted in Appendix [Al where we consider the interac- 
tion of test particles with one ideal wave. 

6.4. Heating rate in test-particle simulations 

In order to make a quantitative comparison between 
our analytical heating rates and our test particle cal- 
culations, we will use the explicitly parameterized dif- 
fusion coefficients of Equations [TU1 Q21 and [121 cor- 
responding respectively to our exponential-decorrelated 
TTD, Gaussian-decorrelated TTD, and FTB models. 
Throughout this section, we will generally assume all 
species are at a constant reference temperature of kT = 
m s Vg. Thus the important scaling is with respect to v s 
(or m s , equivalently) . We choose 7 = 1 and 7 = 0.5 
for the exponential and Gaussian decorrelation models, 
respectively, as these values provide a reasonable fit to 
our simulations. 

We simulate distributions of particles that are 
Maxwellian in v = («j_ + v 2 ) 1 / 2 , and calculate their heat- 
ing rate E = [Ef — Ei)/(tf — U). We normalize the 
initial energy density of the distribution to 3/4 pc 2 , so 
that the test particles represent the energy of the ions 
or electrons in a proton-electron isothermal MHD fluid. 
(Because we use test particles, all of our heating rates 
may be straightforwardly adapted to a different particle 
density by multiplying by a factor ptest/p-) We calcu- 
late a numerical heating rate by fitting a straight line 
to the test particle energy as a function of time, from 
t = 0.25- lL/c s . 

We then determine the coefficients {C\, C2, C3, C4} in 
the analytic models for D\\ by comparison to the test par- 
ticle diffusion coefficients. For each value of e, we focus 
on the value of v± where our Gaussian analytic model 
best matches the location of the high-wi cutoff in the 
parallel diffusion coefficient (in e.g. Figure [71 this cutoff 
is around V\i ~ 2 c s ). Then we choose the to match the 
normalization. The normalization of the analytic curve 
in Figure [Jj is a result of this procedure. 

Finally, we calculate the associated "analytic heating 
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TABLE 2 

Summary of dimensionless constants C, [fi = 1)0 

a These c onstant s n orm alize the diffusion coefficients given by 
Equations IT01IT31 and [19] 





k = 0.01 


e = 0.1 


i = 1 


Ci 


0.085 


0.07 


0.065 




0.15 


0.12 


0.11 


c. 


0.25 


0.17 


0.2 


C A 


0.5 


0.3 


0.4 



rates" by numerically evaluating the integral in Equa- 
tion [20] (integrating the diffusion equation over parallel 
and perpendicular velocities) using our fits for the di- 
mensionless coefficients Ci. We expect the Ci to depend 
only weakly on parameters of the turbulence such as e, 
L, etc. In Table [2] we provide our approximate values 
for these normalization coefficients for three runs at e = 
0.01, 0.1, and 1 c? s /L, all at f3 = 1; we discuss the depen- 
dence on j3 in §6.51 Over a factor of 100 in e, the Ci do 
not change significantly. 

Figure [5] shows how the analytic heating rates com- 
pare to our test particle results for several different val- 
ues of e. Our numerically calculated heating rates never 
asymptote to a constant at high-u s . It is thus clear that 
the Gaussian decorrelation function (orange curves) pro- 
vides a better fit to the simulation data (black curves) 
in the high-t> s regime, where it scales as ln(v s )/v s , as 
opposed to the exponential decorrelation prescription 
(blue curves), which is independent of v s . Similarly, 
the delta function heating rate (pink curves) are typi- 
cally too steep. In the high velocity regime, the heat- 
ing rates are reasonably well-fit by a functional form 
e par t — 0.33 v\{v s /c s ) A , though this expression only 
applies for j3 = 1. 

The dominant heating mechanism in Figure |9] depends 
on the thermal velocity of the test particles. Above the 
linear theory resonance at «|| ~ v s \ ow , the TTD contri- 
bution to the diffusion coefficient scales like /i 2 / 1 ^ <x 
v j_/ v ^ where a — 2 or 3 depending on the decorrelation 
model. Thus, in the bulk of an isotropic thermal distribu- 
tion, where v± ~ uii ~ v s , TTD is increasingly important 
at higher v s . For low thermal velocity particles, on the 
other hand, /i — > and FTB dominates. 

As discussed in f|4] linear theory implies an asymptotic 
heating rate E oc v" 1 for high v s . This in turn implies 
that electrons, with thermal velocity v e = ^mp/m e v p ~ 
43 v p (where m p and v p are respectively the mass and 
thermal velocity of protons) are heated much less effec- 
tively than protons by the compressive fluctuations in 
(3 ~ 1 MHD turbulence. In Tabled we provide estimates 
of E p / (E p + E e ), the proton-to-total heating ratios in our 
test particle simulations, for a range of e and (3. We as- 
sign protons a value of v s — c s /\/2, appropriate for an 
equal-temperature electron-proton plasma. These calcu- 
lations indicate that the electron heating rate is typically 
smaller than the proton heating rate by a factor of 2-5, 
rather than 43. This is primarily due to the asymptotic 
\n(v s )/v s scaling of the simulated heating rates, corre- 
sponding to our Gaussian decorrelation model (see eqn. 
[35]). It is interesting to note that the proton-to-total 




Fig. 9. — Test particle heating rate as a function of test parti- 
cle thermal velocity v s , for 0=1 and several different values of 
the turbulence driving rate e; we also show a comparison to our 
analytically predicted heating rates from i[4] Black curves corre- 
spond to test particle results, blue curv es correspond to an expo- 
nential decorrelation model (see eqn . I24jl , orange curves correspond 
to Gaussian decorrelation (eqn. 1 25 D . and pink curves to a purely 
linear theory diffusion coefficient. In each case, FTB heating dom- 
inates at low v s , while TTD dominates at high v s . The Gaussian 
decorrelation model provides a notably more accurate scaling at 
high v s than the exponential model. At lower e, the delta-function 
heating rate becomes more plausible; this is consistent with the 
idea that the turbulence is weaker in these simulations. Test par- 
ticle heating of e.g. protons with v s = c s /\/2 is wi thin a factor of 
4 of the turbulent cascade rate e in all cases. See t|6.4l for a more 
detailed discussion. 



TABLE 3 

Proton-to-total heating ratio in test particle simulations 



E p /(E P + E e ) e{cj/L) f3 



0.84 


1 


1 


0.84 


0.1 


1 


0.85 


0.01 


1 


0.81 


0.1 


3 


0.74 


0.1 


0.3 


0.69 


0.1 


0.1 



heating rates we find are consistent with empirical infer- 
ences of proton vs. electron heating in the solar wind 
(|Cranmer et al. 2009|), both in magnitude and in the in- 
creasing electron heating for smaller /?. 

In Figure [HI our analytic calculations with Gaussian 
resonance broadening overestimate the magnitude of the 
TTD heating at high v s , particularly for smaller e. This 
is for two reasons. The first is that we fit heating rates 
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Fig. 10. — Test particle heating rate as a function of test particle 
thermal velocity v s for several different values of plasma f} at fixed 
kinetic energy input rate e = 0.1 c 3 s /L. 

over the baseline t — [0.25, 1] L/c s . However, because the 
distributions are driven away from isotropy, the heating 
rate becomes less efficient over time, and so the test par- 
ticle energy increases somewhat sub-linearly. Our an- 
alytics assume instead an isotropic Maxwellian at the 
original thermal velocity. 

A more important effect is that the analytic model used 
in Figure[3]is purely strong turbulence, with ui n \ ~ wunear 
(see £|3.2.1j) . However, the turbulence in the simulations 
with lower e is in fact weaker, with cj n i < wii nC ar- We 
consider the diffusion coefficients resulting from a com- 
bination of weak and strong turbulence in Appendix [Bj 
this results in a sharper resonance which approaches the 
delta- function of linear theory in the Ma — > limit. 
Thus we might expect a heating rate closer to the lin- 
ear theory heating rate (the pink curves in Figure [9]) for 
the runs with smaller e, as indeed is the case in Figure 
[SJ However, this interpretation is complicated by the 
fact that we do not find a clear peak in the test parti- 
cle velocity diffusion coefficients at low e (see Figure HJ. 
This is likely due to phase-decorrelation broadening, as 
discussed in ^3.2.41 

6.5. Dependence on /3 

Figure [10] shows the simulated heating rates for a ther- 
mal distribution of particles for several different values 
of plasma /3 at fixed e = 0.1 c^/L, measured over a longer 
baseline of t = [1, 5] L/c s . We do not attempt to make a 
quantitative comparison with our analytic model. How- 
ever, qualitatively, three effects are clear. First, FTB 
heating at low v s decreases in effectiveness at low /3. This 
is due to the decreasing curvature of the typical magnetic 
field line involved in FTB interactions: vl/va decreases 
as /3 decreases at fixed sound speed and e. 

Second, the contribution of slow modes to TTD de- 
creases with decreasing beta. The fraction of the slow 
mode energy in parallel magnetic field compressions is 
oc /3/(/3 + 1). Thus, at low /3 TTD heating due to the 
slow modes becomes less important. 

Finally, for f3 = 0.1, there is a clear bump in the heat- 
ing rate at high thermal velocities, which moves to lower 
v s at higher fi. We associate this peak with TTD heating 



by fast modes, which have a phase velocity approaching 
va in the low-/3 limit. The reduced heating in the j3 = 0.1 
case relative to the /3 = 0.3 is simply due to the fact that 
our simulations have a greatly reduced proportion of fast 
mode energy at lower /3, by a factor of ^3. 

6.6. Dependence on driving correlation time 

Figure [TTJ shows the heating rate of a thermal dis- 
tribution for different i C orr, the correlation time of the 
Ornstcin-Uhlcnbcck process with which we drive the tur- 
bulence. Over a factor of ^400 in i corr , there is only a 
factor of ~ 4 change in the heating rate; this suggests 
that to zeroth order, our results are insensitive to the 
value of t corr . However, when i corr — > £mhd (not shown 
in Fig. [TTj) . where iMHD is the MHD timestep, we find 
that our results approach the limit of an uncorrelated 
driving scheme (typically, iMHD ~ 8 x 10~ 4 L/c s in our 
highest resolution simulations). In this limit, we find 
acceleration mimicking the cyclotron resonance for high- 
gyrofrequency particles which should not be resonant, 
because of the artifi cially high frequen cies introduced by 
the driving (see also lLehe et al.l[2009f ). 

Increasing i corr from this minimum value while fixing 
other properties of the turbulence systematically affects 
the turbulent kinetic energy in each of the MHD modes. 
The other physically relevant timescale in our calcula- 
tions is the outer-scale eddy turnover time, which for 
our fiducial simulation (see Table [JJ is approximately 
ieddy ~ L/vl, which is approximately given by ~ 3L/c s 
for the fiducial c ase. Again applying the Fourier spectral 
decomposition of lCho fc Lazarianl (|2003f ). increasing i corr 
from 0.01 L/c s to 4.0 L/c s decreases the Alfvenic energy 
from 55% to 45%, decreases the fast mode energy from 
15% to 3%, and increases the slow mode energy from 
30% to 45% (while decreasing the overall kinetic energy 
in the turbulence by roughly 30%). Note that we hold 
the input power e fixed as £ cor r varies. 

We also find that despite driving no modes with paral- 
lel wavelengths longer than L/2, long- wavelength modes 
with k\\VA ~ i^orr naturally appear in the developed 
turbulence. We believe that this is due to frequency- 
matching between the correlation time of the driving and 
the natural frequency of long-parallel-wavelength modes 
(see Appendix |C|) . 

These varying proportions of energy affect the thermal 
heating rate in ways which are largely consistent with our 
interpretation of the heating (see Fig. [TTj) . The decreas- 
ing energy in fast modes manifests itself in a factor of ~ 4 
decrease in the TTD heating rate at high thermal veloc- 
ity, while the increase in slow mode energy contributes 
to increasing the heating rate just below v s ~ c s . The 
decrease in FTB heating is consistent with the decrease 
of vl at longer correlation times. 

6.7. Non-thermal acceleration? 

In previous sections, we have focused on short- 
timescale interactions between particles and turbulence 
to calculate diffusion coefficients and corresponding heat- 
ing rates. On longer timescales, it is unclear whether 
these interactions produce a thermal or non-thermal evo- 
lution of the distribution function. To investigate this 
question, we carried out test particle simulations lasting 
for longer times. 



14 




ltr 1 io° 10 1 

v a (c s ) 

Fig. 11. — Test particle heating rate as a function of test particle 
thermal velocity v s for several different values of the OU turbulent 
driving correlation time, t CO rr, in units of L/c s . All other parame- 
ters are our fiducial parameters, listed in Table [l] Changes in the 
heating rate are largel y du e to varying proportions of energy in 
the MHD modes (see q6.7l) . but the overall results are consistent 
at the factor of ~ 4 level over a factor of ~ 400 change in t CO rr ■ The 
large-scale eddy turnover time corresponds to t CO rr ~ a few L/c s . 
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Fig. 12. — The distribution function of an initially thermal dis- 
tribution of test particles with v s = c s evolving in turbulence with 
ft = 1 and e = 0.1, plotted at t = (solid curve) and t = 20 L/c s 
(dotted curve); this corresponds to ~ 10 eddy turnover times of 
the turbulence. We also plot a thermal distribution with the same 
energy as the final distribution (long dashed curve). The distri- 
bution becomes mildly non-thermal due to the turbulent diffusion, 
but TTD does not lead to a power-law tail. 

Figure [T2] shows the distribution function (dotted 
curve) resulting from a simulation with e = 0.1 (? s jL run 
for t = 20 L/c s , corresponding to several eddy turnover 
times. The thermal energy of the particles has increased 
by ~ 50%. For comparison, Figure [12] also shows a 
Maxwellian with the same energy as the final distribu- 
tion (long dashed curve). The final distribution is clearly 
non-thermal, in that there is somewhat more energy in 
high- velocity particles, and less for low-velocity particles. 
However, the distribution is only weakly non-thermal, in 
that there is no evidence for the formation of a power-law 
tail at high velocities. 



This weak non-thermality is not that surprising, as 
even our broadened TTD resonance diffusion coefficient 
is largest for particles with vu ~ va- Our result is in- 
consistent with recent work claiming that compressible 
turbulence generically leads to non-thermal power law 
tails f(v) oc ^~ 5 , as observ ed in the solar wind (see 
iFisk fc Gloecklerll2007ll2008t ). However, their formalism 
assumes an isotropic or nearly isotropic distribution func- 
tion, enforced by particle collision s or microscale insta - 
bilit ies, and has been critici zed by IJokipii fc Led (|2010T > 
and ISchwadron et al.l (|2010f) . among others. We will in- 
vestigate this point in more detail in future work. 

7. CONCLUSIONS 

We have studied the interaction between charged test 
particles and low-frequency, large-scale MHD turbulence. 
This interaction is important in a wide range of astro- 
physical systems, including the solar wind and cosmic ray 
transport through the galaxy. The coupling between par- 
ticles and MHD turbulence leads to velocity-space evo- 
lution of the particles, including diffusion, heating, and 
acceleration. We have used simple physical arguments to 
motivate analytic models of resonance broadening due 
to the interaction between particles and strong turbu- 
lence f £|3.2.1[) . Furthermore, we have shown that the 
non-resonant interaction of charged particles with mov- 
ing, curved magnetic field lines (see Fig. [IJ is important 
for a full understanding of the velocity-space diffusion of 
particles in a turbulent plasma f £)3.2.3[) . 

We have calibrated these analytic models of velocity 
space diffusion against simulations of charged test par- 
ticles in fully dynamical MHD turbulence (Q. These 
calibrations are summarized in Equations I10l and ll31 and 
Table [5J We anticipate that these calibrations of the ve- 
locity space diffusion of test particles in MHD turbulence 
will be useful for a wide range of future astrophysical and 
hcliospheric applications. 

Our most important results include: 

• The transit-time damping resonance is highly 
broadened in MHD turbulence, relative to the delta- 
function prediction of linear theory. Our phe- 
nomcnological model for resonance broadening, 
which describes wave decoherence in strong MHD 
turbulence, generically leads to a velocity diffusion 
coefficient which approaches a constant at low-vy 
and a high-«|i power-law tail (see Fig. [2]). We also 
see evidence for phase-decorrelation broadening, in 
which finite-amplitude turbulence accelerates par- 
ticles into and out of resonance over relatively short 
timescales (see H3.2.4[) . The significant broadening 
we find implies that many particles, not just par- 
ticles moving with the phase velocity of the waves, 
can strongly interact with the turbulence. Presum- 
ably this conclusion would also apply to the Lan- 
dau resonance in turbulence with parallel electric 
fields, although in this study we have limited our 
considerations and simulations to ideal MHD. 

• Heating rates for high thermal velocity particles are 
inconsistent with a slow mode decorrelation model 
which is exponential in time (i.e., a Lorenztian res- 
onance function). A Gaussian model for slow mode 
decorrelation produces heating rates at high ther- 
mal velocities which are much more consistent with 
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our test particle calculations (see Fig. . Because 
of their weaker nonlinearity and longer decorrela- 
tion time, we are not able to distinguish between 
the exponential and Gaussian decorrelation models 
for fast modes. 

• Fermi Type-B interactions, wherein particles are 
slung around by moving, curved magnetic field 
lines (see Fig. [1} , are critical in describing the full 
velocity space diffusion of charged particles with 
MHD turbulence (see g3jO] and Fig. B). These 
interactions have a non-resonant character, and ac- 
celerate particles independent of the particles' mag- 
netic moment fj,. In general, FTB dominates for 
particles with low fi (where the TTD interactions 
become correspondingly weak) or for particles with 

1/2 

vn ^> (va/vl) v± (such that even the broadened 
TTD resonance has fallen off significantly) . 

• FTB dominates the heating of particles with ther- 
mal velocities v s much less than the fast and slow 
mode phase speeds, while TTD dominates the 
heating of high-u s particles (see Fig.[S]). FTB heat- 
ing is thus particularly important for minor ions 
which have thermal velocities less than the plas- 
mas sound speed. For our fiducial case, with Ma 
similar to the solar wind, FTB and TTD contribute 
a similar amount of heating to protons. 

• TTD can efficiently damp the turbulent energy in 
compressive MHD modes: we find test particle 
heating rates comparable to the turbulent energy 
cascade rate for a wide range of plasma param- 
eters (see Fig. |9]). Furthermore, electron heating 
is comparable to proton heating for the range of 
/? ~ 0.1 — 3 we studied (see Tabled]). Our estimated 
proton-to-total heating ratios are consistent with 
empi rical studies in the solar wind (jCranmer et al.1 
2009). We note, however, that our calculations do 
not include processes that damp the Alfvenic com- 
ponent of the turbulent fluctuations at small scales 
and so do not capture all of the heating that is 
likely important in the solar wind and other astro- 
physical plasmas. 

• MHD turbulence does not efficiently accelerate col- 
lisionless test particles out of the bulk of a thermal 
distribution (see Fig. [T2"]) . We find no evidence for 
the formation of a power-law tail even after many 
turnover times of the turbulence on large scales. 
Instead, most of the turbulent energy is converted 
into thermal energy of the bulk of the plasma. 

In our simulations, ~ 45 — 50% of kinetic energy is in 
slow modes, while ^3 — 5% is in fast modes. As men- 
tioned in fJSJ this is notably higher than in the near-Earth 
solar wind, where ~ 10% of energy i s in slow modes and a 
negligible fraction is in fast modes (|Howes 
Because compressive wave modes heat test particles so 
efficiently (at a rate generally comparable to the turbu- 
lent cascade rate e), any energy initially in those modes 
would be quickly damped out. This is consistent with the 
fact that the solar wind contains a smaller proportion of 



compressive energy in the inertial range than naive ideal 
MHD simulations. 

We also find that the Alfvenic component of the turbu- 
lent cascade will be significantly damped by Fermi Type- 
B interactions at the outer scale where the turbulence 
spectrum is nearly isotropic. However, FTB interactions 
quickly become weak on smaller scales, and so their effect 
on the the inertial range of the turbulence can probably 
be neglected. Transit-time damping, on the other hand, 
damps energy out of all decades in wavenumber at equal 
rates, but requires energy in compressive modes to be 
effective. 

Another interesting consequence of FTB interactions 
is that velocity diffusion depends only on the turbulence 
properties, and is independent of particle mass or charge. 
This implies that the heating time, i/j ~ n^T/fi, is in- 
dependent of particle mass and charge. The cooling or 
expansion time t c ~ r/v sw in the solar wind is also in- 
dependent of particle mass and charge (where r is the 
heliocentric distance and v sw is the expansion velocity 
of the solar wind). If i/j 3> t c , the particles will sim- 
ply cool by adiabatic expansion, and their final tem- 
peratures will be determined by their initial tempera- 
tures. However, if th -C t c , then the particles will quickly 
heat up until th ~ t c . If we assume that FTB interac- 
tions are the most important heating process for minor 
ions, this balance implies a temperature which is inde- 
pendent of charge but proportional to the species mass, 
T s ~ m s vlr/(k B Lv A v sw ). 

The interaction between particles and compressive 
MHD turbulence has been invoked to explain the v ~ 5 
distribution at high velocities observed in the solar wind 
(jFisk fc Gloecklerl 120071 . 120081 ). Our results are not con- 
sistent with these models in that we see no evidence for 
the development of a power-law tail to the distribution 
function even after many turnover times on large scales 
(see Fig. [12] and g67l). 

In future work, we intend to implement simple pitch- 
angle scattering of test particles, mimicking the effects 
of small-scale plasma instabilities such as the firehose 
and mirror instabilities. It will be interesting to assess 
whether our results on particle acceleration and the long- 
term evolution of the distribution function change in the 
presence of significant pitch-angle scattering. In addi- 
tion, our test particle methods are sufficiently general 
that they may be applied in simulations of turbulence 
which are relevant on smaller scales, e.g. Hall MHD tur- 
bulence, to probe the gyroscale transition in the solar 
wind. 
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rial is based on work supported by the National Sci- 
ence Foundation Graduate Research Fellowship under 
Grant No. DGE- 11 06400. Additionally, this work was 
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0752503. Computing time was provided by the National 
Science Foundation TeraGrid/XSEDE resource on the 
Kraken and Frost supercomputers. Additional computa- 
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APPENDIX 

A. LANDAU- TYPE RESONANCE BETWEEN TEST PARTICLES AND ONE WAVE 

To consider the effects of particle trapping and bin-crossing on our method for measuring velocity diffusion coef- 
ficients, we have simulated a simpler self-contained one-dimensional toy problem, where particles feel a sinusoidal, 
travelling acceleration of the form a = ciq sin (kx — kv p t). This is effectively the Landau problem with test particles. 
(The toy problem is independent of Athena, but we have confirmed that the results are qualitatively identical to test 
particles interacting with a single slow mode wave in Athena.) 

The mean-square change in the velocity of test particles {(5v) 2 ) for several waves with different amplitudes is shown 
in Figure [13] We plot ((Sv) 2 } because this quantity converges to a constant at long times. If each bin in v were truly 
a delta- function, and there was no numerical loss of accuracy, then {(Sv) 2 } for test particles in a given bin would reach 
a maximum after approximately 0.5L/(v — v p ) (corresponding to half of the wave-particle interaction time), and then 
decrease back to zero as every particle returned to its original phase with respect to the wave, though offset by 2ir. 
However, because particles are smoothly distributed in velocity within each bin, slight initial velocity differences lead 
to the destruction of this phase coherence over time, and ((Sv) 2 ) converges to a constant value of roughly half its 
maximum. 
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Fig. 13. — Mean-square velocity change for ID toy problem of test particles interacting with a wave of phase speed v p , measured after 
60 wave periods (long enough that the plotted quantity has converged). The ratio of the linear wave time to the nonlinear acceleration 
time (evaluated for particles at the resonance) ranges logarithmically from 0.001 (solid curve) to 1 (dashed curve). For small-amplitude 
waves, a sharp resonance at v p is evident. As the amplitude grows larger, particle trapping effects associated with finite amplitude lead 
to a broadening of the resonance. A similar effect occurs in MHD tur bulenc e, where finite-amplitude fi\7B forces move resonant particles 
off-resonance. We call this effect phase-decorrelation broadening (see i|3. 2.41 1. 



Acceleration by a small-amplitude wave produces a very sharp resonance around the wave phase velocity, while 
particle trapping significantly broadens the resonance for large-amplitude waves. For this problem, a small amplitude 
wave is one for which tu a /t acc -C 1, where t\i n ~ l/v p k is the linear wave period and t acc ~ v/o,q is the time for a 
particle's velocity to change significantly due to the wave-particle interaction. The transport in velocity space is not 
truly diffusive for this toy problem, except for particles very near the resonance. 

The results of this simple problem are directly analogous to test particles interacting with MHD turbulence. When 
iiin/iph ^ 1 (where i pri is the time required for particles to significantly move out of phase with a wave), particles 
experience significant acceleration and quickly fall out of phase with a previously resonant wave, leading to strong 
resonance broadening. The situation in turbulence is discussed in ^3.2.41 

B. RESONANCE BROADENING IN WEAK TURBULENCE 
B.l. Turbulence Model 

In the main text, we considered a resonance broadening model where the underlying turbulence spectrum was purely 
strong GS turbulence. In this appendix, we discuss an expansion of the model to the case of weak turbulence. 

We consider MHD turbulence that is driven isotropically with an outer scale eddy velocity vl < va- The turbulence 
at the outer scale is weak in the sense that the linear time scale for wave propagation is shorter than the nonlinear 
time scale for eddy turnover, L\\/va < L±/vl, where isotropy implies that Ln = L± = L. This implies that waves live 
for more than one wave period before they decay. We may express this weak turbulence requirement more generally 
by noting that the non-linearity parameter x = ViI\\/vaI±. < 1, where l\\ and l± refer to scales smaller than L, and vi 
is the turbulent eddy velocity on that scale. 

However, weak turbulence cascades only in the perpendicular direction to higher k± (Sheb alin et al.l 119831 : 
iGaltier et al.l [2000), while no structure on smaller parallel scales develops. This implies a ID power spectrum 

Ek oc fcJ 2 (|Ng fc Bhattacharjed Il997h . However, x ^ k± 2 , so that the non-linearity increases at smaller scales. 
Eventually, when % ~ 1, the turbulence reaches the strong critically balanced state of GS, and the cascade proceeds 
to smaller scales at fixed x 01 order unity. Critical balance refers to the balance of wave and eddy time scales. This 



implies that most of the turbulent energy is contained in eddies with k 



'2/3 



and a ID power spectrum Ek oc k 



-5/3 



The transition between the weak and strong regimes occurs at a wavenumber k± t , 
is the Alfvenic Mach number at the outer scale. 
More quantitatively, we define the 3D power spectrum as 

7(k) = J w (k)+/ S (k), 

where the weak and strong contributions to the power spectrum are given respectively by 



(2%/L)Mj 2 , where M A = v L /v A 



I w (k) = A^r^u - 2n/L)g 



and 



J.(k) = A s kl 10/3 g 



k\\L 
~2T 



k±, c 



1-5 



(Bl) 



(B2) 



(B3) 
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where A w and A s are dimensional normalizing coefficients defined below, and g(x) will be treated as a step function, 
equal to 1 if \x\ < 1, and otherwise. g(x) is used in two ways: the first g(x) in Equation IB3I accounts for the fact that 
power only resides inside the critically balanced cone of wavenumbers, while the other instances model the transition 
from weak to strong turbulence at k± — fcj_ iC . 

We fix the normalizing coefficients A w and A s by first requiring that the total energy in the turbulence be given by 
an integral over the power spectrum, 

' ^d 3 k/(k). (B4) 



1L 
2 



(2tt) 3 



Because there are two coefficients A w and A s , we require another condition, which is given by the requirement that 
the ID weak and strong power spectra must match up at the transition scale: 



dk\\I w (k) 

Taken together, these two requirements imply that 



fc_L,c 



dk\\I a (k) 



(B5) 
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where F is given by 



A s = (2irf^M A 
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and N = L/L m [ n is the range of scales in the turbulent cascade. For N ^> 1, F is approximately constant and 
equal to 1 until Ma is near unity, and then increases sharply to 3/2 as Ma —> 1- Recent work which demonstrates a 
transition from weak to strong turbulence (jHowes et alJl20lTbt iVerdini fc Grapp"m1l2012| ) supports the basic form of 
this turbulence model. 



B.2. Resonance Broadening 

We use the same resonance broadening model as in the main text (see Equation [5]) . For weak turbulence, this 
approach must be modified slightly to reflect the fact that while \ is the fractional energy change of a wave in one 
wave period, the energy change could be randomly positive or negative, and the energy of the wave will diffuse in 
energy space. Thus we estimate that the actual decoherence time of waves will be typically given by Wdcc = X 2 ^iin- 

As before, the exponential resonance function is given by 

7 2 ^dec + fc |( V ll ±V p) 2 

where Wdcc scales differently with k± in the weak and strong cases. We see that the resonance function is still peaked 
at |f||| = v p , but becomes a Lorentzian in vu rather than a delta-function. 

The diffusion coefficient is a sum of the contributions from the weak and strong turb ulenc e components. To calculate 
the contribution from weak turbulence, we use the weak power spectrum of Equation IB21 and a non-linear frequency 
given by Wdcc = k^v 2 L /vA to find 

G [u\ ( M\+u 2 \ ( Mi , 
D\\ <w = T^-T -4 M 771773 ^ 77TT-, (BIO) 



%w - 16ul \ ul { M\ (u 2 _ +1)1 \MX( 




[a 




where 



a . i^^Mi, (bh, 



and u± = {v\\ ± v p )/(^v A ). 

Beyond the transition to strong turbulence at fc_i_ iC , the decorrelation frequency is given by Wdcc = = 
{2ttva/ L)(k±/ k±, c ) 2 / 3 , and the power spectrum is given by Equation IB3I In this case we find 



D \\,s = ln [ nm a] I f 1 + ~f) — — (~r arctanu_ + arctanu + j \ . (B12) 
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which is again strongly broadened with an asymptotic form as 3> v p of D^ s cx v>, 2 . For Ma < 0.5, the strong 

turbulence diffusion c oeffici ent in Equation IB12I is significantly more broadened than the weak turbulence diffusion 
coefficient in Equation IB101 as we show explicitly in Figure [T4l discussed below. 

As in the main text, we also examine a Gaussian decorrelation model, where we replace the decorrelation term 
7Wdec^ with (7Wdec^) 2 - I n this case, we find the Gaussian resonance function, 



R 



gauss 



(k) 



27<^dc 



■ exp 



k \\ v \\) 7 

4 7 2 ^dcc 



so that the (^-function becomes a Gaussian resonance. Now, using Equation [31 we find 



(B13) 



and 
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where erf (a;) is the error function. In the 3> v p limit, -Dii^ cuts off exponentially, while D^ s 




Fig. 14. — Diffusion coefficients (arbitrary normalization) for transit-time damping by slow modes in exponential (left) and Gaussian 
(right) decorrelation models. Ma = Vl/va = B^/Bq = 0.25, and v p ~ 0.81, Va = v2, corresponding to slow modes with /3 = 1. For this 
choice of Ma we see that weak turbulence on the outer scales is the dominant contribution to the velocity diffusion, although the |f||| — 3 
tail from strong turbulence is important at high jjii in the Gaussian model. 



In Figure [T4l we plot representative examples of Du for exponential and Gaussian decorrelation models, with the 
same arbitrary normalization. We choose 7 = 1 and 7 = 0.5 for the exponential and Gaussian models respectively, 
as in the main text. The diffusion coefficients display a sharp peak around the resonant phase velocity. In Figure [5] 
in the main text we plotted the corresponding diffusion coefficients for a run with even lower amplitude turbulence 
{Ma — 0.11), for which this model would predict a yet sharper resonance. No clear peak could be seen. Wc interpret 
that this is the result of additional phase-decorrelation broadening due to finite-amplitude turbulence ( £I3.2.4[) . which 
we do not quantitatively model. 

C. TURBULENCE POWER SPECTRA 

In Figure 1151 we show the 2D power spectra of fully saturated turbulence for two different simulations, one with 
^corr = 0.1 L/c s (left panel) and one with i corr = 1.5 L/c s right panel. In all other respects the simulations are identical 
and have the fiducial parameters of Table [1] with an eddy turnover time of approximately i ddy ~ (L/2)/vl ~ 1.5L/c s . 

The simulation with the longer correlation time develops significant power on low fcii , which are not driven by our 
turbulence driving routine, while the simulation with the shorter correlation time does not develop power on similar 
scales. We believe that this is due to frequency matching of waves with the frequency content of the Ornstein-Uhlenbeck 
driving routine, so that modes with k\\VA ~ naturally develop power despite not being actively driven in k-space. 

We also note that the longer-£ corr simulation has a one-dimensional power spectrum consistent with P cx /c~ 5 / 3 , while 
the simulation with the shorter correlation time has a somewhat steeper power spectrum, possibly consistent with 
P cx fc -2 (though our lack of inertial range does not allow us to make this statement with any confidence). This 
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k\\L/47T k\\L/47T 

Fig. 15. — Two-dimensional plots of P(k)fc 3 in {fcy , fcx}-space for saturated turbulence with t corv = 0.1 L/c s (left panel) and t CO rr = 
1.5L/c s (right). All other parameters are for our fiducial case, given in Table [T] Color normalization is logarithmic and identical, with 
white and yellow corresponding to the highest values, and red and black corresponding to the lowest values. Despite the fact that modes 
with fey < Att/L are not driven, frequency-matching of long- wavelength Alfven waves with the long t C orr leads to the development of power 
in low-fc|| modes. In the run with a shorter correlation time, power is much more concentrated in the driven volume of fc-space. See 
discussion in Appendix [Cl 

seems consistent with the idea that the extra power at low k\\ in the run with the longer correlation time provides 
anisotropy that can actually lead to critical balance on the outer scales. In other words, \ ~ k±VL/k\\VA ~ 1 because 
ku /k± ~ vl /va as a result of frequency matching. We do not attempt to explicitly model this behavior in the main body 
of the paper, however, choosing instead to focus on the simpler and more generic case of isotropic strong turbulence. 
Furthermore, the test particle diffusion coefficients do not demonstrate any obvious peakiness for the shorter t COII 
case, suggesting that despite the apparent weaker turbulence, the diffusion coefficients are not well-described by the 
turbulence model we discuss in Appendix [B] We believe that this is a result of the phase-decorrelation broadening 
discussed in i)3.2.4l 



